Geometric scaling of purely-elastic flow instabilities 
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We present a combined experimental, numerical and theoretical investigation of the geometric 
scaling of the onset of a purely-elastic flow instability in a serpentine channel. Good qualitative 
agreement is obtained between experiments, using dilute solutions of flexible polymers in microflu- 
idic devices, and two-dimensional numerical simulations using the UCM model. The results are 
confirmed by a simple theoretical analysis, based on the dimensionless criterion proposed by Pakdel- 
McKinley [l|, Q] for onset of a purely-elastic instability. 



PACS numbers: 47.50.-d, 47.61.-k, 47.20.Gv 

Purely-elastic flow instabilities occur for low Reynolds 
number flows of viscoelastic fluids. Such elastically- 
driven instabilities - which are entirely absent for the 
equivalent Newtonian fluid flow - have been experimen- 
tally observed in a range of macroscopic flow geometries 
such as viscometric Couette or plate-plate devices 0, H| . 
The onset of instability has been attributed to a bal- 
ance between normal stresses and streamline curvature 
0. Pakdel and McKinley (Pak-McK) Q, proposed 
an elegantly simple dimensionless criterion to unify the 
experimental observations up to that date. However, be- 
yond the lid-driven cavity results used in [l[ , systematic 
studies of the instability threshold as a function of the 
curvature of a given flow are still lacking. In this Let- 
ter we investigate the geometric scaling of the onset of 
a purely-elastic instability in serpentine channels of rect- 
angular cross-section, combining experiments, numerical 
simulations and a simple theoretical analysis. In such 
wavy channels the geometric influence on the curvature 
is dominant and the curvature is essentially independent 
of flow rate. A major complication of investigating cur- 
vature effects in most other complex flows - e.g. abrupt 
contractions or lid-driven cavity flows - is that the dom- 
inant flow curvature is not solely determined geometri- 
cally but by the precise dynamics of the flow itself (e.g. 
varying locally and with flow rate). These elastic insta- 
bilities have proven to be a good way to obtain efficient 
mixing in microfluidic devices Newtonian fluid flow 
in microfluidic devices occurs typically in the laminar 
regime, due to the characteristic small length scales, and 
mixing occurs primarily via diffusion, a very slow pro- 
cess. However, the addition of small amounts of flexible 
polymers with high molecular weight is sufficient to im- 
part strong non-Newtonian behavior, and can even lead 
to so-called "elastic turbulence" [f§, with increased flow 
resistance and enhanced mixing performance Q. It is 
thus of great importance to understand the underlying 
mechanism that leads to the onset of instability and test 
the robustness of the Pak-McK criterion. 




FIG. 1: Sketch of the microchannel; typical instantaneous flow 
patterns in a dilute solution of PEO. From top to bottom: 
stable flow; slightly unstable flow (close to onset of elastic 
instability); unstable flow. 



We use microfluidic devices, which are important in 
applications such as lab-on-a-chip, fabricated in poly- 
dimethylsiloxane (PDMS) from an SU-8 mould, us- 
ing standard soft-lithography techniques, which allow 
straightforward variation of the flow geometry. Fur- 
thermore, due to the small dimensions of the channel, 
high shear rates are easily achievable, thus enhancing the 
non-Newtonian behaviour observed, while keeping iner- 
tia small. The experiments were performed in serpentine 
channels with varying radii of curvature and channel di- 
mensions, using dilute solutions of a flexible polymer in 
different solvent viscosities. The onset of a purely-elastic 
instability (see Fig. [1} was investigated function 
of the relaxation time of the flexible polymer solution 
and the geometrical characteristics of the channel (ra- 
dius and width). Two-dimensional (2D) numerical simu- 
lations arc performed using the upper-convected Maxwell 
(UCM) model in identical geometries, assuming that the 
curvature in the 2D base-flow is responsible for the insta- 
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bility (preliminary 3D calculations confirmed the appro- 
priateness of this approach). The combination of both 
approaches, together with a simple theoretical analysis, 
based on the Pak-McK criterion, allowed us for the first 
time to predict the geometric scaling of the onset of a 
purely-elastic instability in curved channel flow. In addi- 
tion our results independently confirm the strength of the 
Pak-McK criterion and open the possibility for extending 
similar analyses to a range of other flows. 

The geometry of the serpentine channels used in this 
study are represented in Fig. HJa) (simulations) and in 
Fig. [T] (experiments). The flow takes place through a 
given number of half loops in a channel of width W (and 
height H) and inner radius R. In the straight parts 
of the inlet and outlet the velocity profile can be esti- 
mated using the analytical solution for fully-developed 
flow in a rectangular channel, valid for fluids with a 
constant shear viscosity (either Newtonian or viscoelas- 
tic). In the curved parts of the channel the flow becomes 
weakly asymmetric and the maximum velocity location 
occurs slightly closer to the inner wall. In between the 
two half loops the flow regains symmetry, before becom- 
ing asymmetric towards the other side in the next half 
loop. The asymmetry becomes less pronounced with in- 
creasing dimcnsionlcss radius R/W. The flow is mainly 
shear dominated, as can be seen on the contour plot of 
the flow-type parameter shown in Fig. [2J} (from simu- 
lations), with limited regions of elongational and rota- 
tional flows close to the centreline where the deforma- 
tion rates are small. The flow-type parameter is defined 
as £ = d + n ' wnere |D| and represent the mag- 
nitudes of the rate of deformation and vorticity tensors 
[§], and varies from £ = — 1 (solid- like rotation), up to 
£ = 1 (extensional flow). A shear-dominated region is 
identified where £ = 0. We define the average shear rate 
as 7 = U/W, with U representing the average velocity 
in the channel. The Reynolds number is here defined 
as Re = pUW/rj, with p and rj representing the density 
and shear viscosity of the fluid, respectively. The Weis- 
senberg number is defined as Wi = A7, where A is the 
relaxation time of the fluid. While Re quantifies the rela- 
tive importance of inertial over viscous forces, the degree 
of elasticity of the flow is quantified using Wi (ratio of 
clastic to viscous stresses). 

The experiments were performed in serpentine mi- 
crochanncls composed by eight half loops (Fig. QJa)). 
We work with channels of different width, W = 60 and 
100 pm, but keep the aspect ratio a = W/H nearly con- 
stant and close to 1.4, a typical aspect ratio for microflu- 
idic devices. The radius of curvature was varied from 
i?=25 to 1950 pm, in such a way that for each chan- 
nel width the dimcnsionlcss radius R/W varies approxi- 
mately from 0.5 to 20. 

Solutions of the flexible polymer PEO (Sigma Aldrich) 
with a molecular weig ht, M w = 2 x 10 6 g/mol, at a con- 
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FIG. 2: a) Geometry and definitions (note: the y-direction 
points to the center along any radial line through each half 
loop); b) flow-type parameter. 



ccntration of 125 ppm (w/w) were used in watcr/glyccrol 
mixtures. The overlap concentration is c* ~ 860 ppm Q 
and we thus work in the dilute limit. The solvent vis- 
cosity T] s varies from r] s = 0.93 mPa.s to 7.4 mPa.s for 
varying concentrations of glycerol at 23°C. For these so- 
lutions the polymer viscosity rj p is approximatively 13% 
of the solvent viscosity, rj s , as measured with a Ubbelohdc 
viscometer. We have also measured the flow curves of the 
polymer solutions (not shown) and no shear thinning is 
observed. For such small concentrations it is very difficult 
to measure the relaxation time directly flp| , therefore we 
estimated A from the Zimm relaxation time. Since the 
flow in the serpentine channel is shear dominated, the 
use of a relaxation time linked to shear is the most ap- 
propriate. From Rodd et al. Q we have calculated the 
Zimm relaxation time for PEO of M w = 2x 10 6 g/mol 
at 23°C in water to be A = 0.36 ms. We have fitted the 
dependence of the Zimm relaxation time on the solvent 
viscosity rj s as found by Rodd et al. [§] (see also 11 1) 
using a power law equation, A = OAr/® 7 , with A given in 
ms and r\ in mPa.s, and use this equation to estimate the 
relaxation times of all the solutions used in the experi- 
ments. 

The experimental protocol can be summarized as fol- 
lows. The solutions are fed to the channel via two inlets, 
and one of them is fluorescently dyed. In this way we 
can visually assess the stability of the flow. The flow 
rate (Q) was varied from 1 to 50 pl/min, and is imposed 
via a syringe pump (PHD 2000, Harvard apparatus). The 
flow is visualized using an inverted microscope (Axio Ob- 
server, Zeiss) coupled to a CCD camera. Starting with 
the lowest flow rate, Q is then gradually increased. Af- 
ter each change in Q we wait a sufficiently long time to 
achieve steady-state flow conditions (on the order of 10 
minutes). The onset of fluctuations in the flow is mon- 
itored at the last loop of the serpentine channel. This 
critical condition defines the onset of the time-dependent 
clastic instability, and the critical flow rate Q c is deter- 
mined. For each experimental condition two separate 
runs are performed using freshly prepared polymer so- 
lutions. We have also visualized the onset of instability 
using fluorescent micro-particles (see Fig. [TJb)). 

A systematic study was undertaken to determine the 
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FIG. 3: Threshold of instability in a channel of width W = 
100 pm and H = 80 jim for a solution of 125 ppm of PEO 
(M w = 2 x 10 6 g/mol) with different solvent viscosities 7] s 
(1.8 mPas (•), 2.8 mPas (■), 4.5 mPas ( ), 7.4 mPas (•«)) 
and W = 60 pm and H = 40 pm (0% of glycerol (-)) (a) 
Critical shear rate -y c as a function of the radius R. (b) Critical 
Weissenberg number Wi c as a function of R/W. 



critical flow rate for onset of an elastic instability as a 
function of the radius of curvature of the channel, i?, 
the width of the channel, W , and the relaxation time 
of the polymer solutions, A. The critical shear rate, 
7c = Qc/{HW 2 ), was determined as a function of R, for 
a first set of channels of width W = 100 pm and height 
H = 80 pm, using different solvent viscosities. The re- 
sults are depicted on Fig. [3] (a) from which it is clear 
that 7 C increases with R. Additionally, one can observe 
that the critical shear rate is higher for lower solvent vis- 
cosities. The critical Weissenberg number, Wi c = A7 C , 
is represented in Fig. [3] (b) as a function of R/W. A 
reasonably good scaling of the data is obtained for all 
sets of experiments. Additionally, the data from another 
set of channels iW = 60 pm, H = 40 pm) are also de- 
picted on Fig. [21(b) for selected solvent viscosities. Good 
agreement between the different channel sizes is obtained, 
therefore confirming that Wi c is a function of R/W. The 
Reynolds numbers for the various experiments presented 
on Fig. [3jvary from 0.2 to 5, therefore, although small, 
inertial effects are not negligible. However, we do not 
observe a modification of the Wi c threshold as a func- 
tion of Re (largest for the small solvent viscosities). For 
the smallest solvent viscosity and a channel of radius 



R = 1950 pm, the flow does not become unstable in the 
range of flow rates accessible in the experiments, and we 
attribute this behaviour to a stabilizing effect due to in- 
ertia for Re > 5 (as also observed in [12|). 

For simplicity, the numerical simulations assume 2D 
planar flow of an incompressible viscoelastic fluid de- 
scribed by the UCM model. The equations that need 
to be solved are those of mass conservation, V • u = 0, 
and the momentum equation, — Vp + V-t = 0, assuming 
creeping-flow conditions (i.e. the inertial terms are ex- 
actly zero). The evolution equation for the extra-stress 
tensor, t, is given by t + ^ T (i) = '77 > where Tm rep- 
resents the upper-convected derivative of r. Despite the 
well-known shortcomings of this simplified model, such as 
the unbounded nature of the steady extensional stresses 
above a critical strain rate (Ae = 0.5), it is the simplest 
differential constitutive equation which can capture qual- 
itatively many features of highly-elastic flows. In addi- 
tion it forms the backbone of many more complex mod- 
els (e.g. the FENE-P, Giesekus and Phan-Thien-Tanner 
models), which simplify to the UCM in certain parameter 
limits, and thus it is extremely general. The governing 
equations are solved using a fully implicit finite-volume 
numerical method, based on the logarithm transforma- 
tion of the conformation tensor, as described in detail in 
[l3| and references therein. 
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FIG. 4: Comparison of experimental and numerical results, 
with the scaling analysis based on Pak-McK criterion. The 
experimental data are the average over the results presented 
in Fig. [3j (b) for the channel of width W = 100 pm. 

Fig. |4]displays the average Wi c , and the corresponding 
standard deviation, of the experimental data, together 
with the predictions from the numerical simulations, as 
a function of R/W. Wi c for the simulations reflects the 
highest Wi for which a steady solution could be obtained. 
In the experiments the ratio r\ v /r\ s is constant and we can 
thus expect qualitative agreement between the experi- 
mental data and the results from the simulations using 
the UCM model. This is indeed the case and good qual- 
itative agreement between experiments and simulations 
is obtained. 
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FIG. 5: a) Sketch of the modified channel; b) streak photog- 
raphy for a dilute solution of PEO: the flow is unstable in the 
part of the channel with small radius (R = 50 um), but stable 
in the part with the large radius (R = 450 um). 



Pakdel and McKinley i] showed that the curva- 
ture of the flow and the tensile stress acting along the 
streamlines could be combined in the form of a dimen- 
sionless criterion that must be exceeded for the onset of 



Til Xu 

rj-y TZ 



> 



Mi nt with n, 



purely elastic instability: 

and 7 representing the local streamline curvature, veloc- 
ity and shear rate, respectively. Tn represents the local 
streamwise normal stress and 777 the local shear stress. 
For the UCM model the normal stress can be written as 
Tn = 2i]X A f 2 . Assuming fully-developed flow conditions 
everywhere (below Wi c ), we can estimate the velocity 
profile assuming a Poiseuille flow, u = §{7[1 — (p^) 2 ], 
with the y-axis centred on the channel and directed along 
any radial line through each half loop towards the centre. 
The shear rate is given as 7 = 12U\y\/W 2 , therefore one 
obtains the following criterion for a purely-elastic insta- 

bility, 36(^) 2 F(y) > M 2 rit , where F(y) = ^~^ ] . 

For each value of R/W we can determine the location 
where the onset of the instability takes place, which oc- 
curs at y c /W where F(y c ) is maximum. Therefore, for 
each value of R/W we can also determine the correspond- 
ing Wi c for onset of instability. In general this has to be 
done numerically, but for the limits of very low and very 
high values of R/W this can be done analytically, leading 
to the following asymptotes and critical locations: 
(i) Wi c -> y c /W \ for § -+ 



'&;y c /W-> 0.289 for § 
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(ii) Wi c -> 0.38M crjt y w , yc/ „ , ^ w 

These asymptotes can be combined in the general form, 

0.38M cri * 



Wi c = 0.38M crit ^0.1+# > which agrees well in the 
whole range of R/W with the exact results that are ob- 
tained by solving numerically the general equation de- 
rived from the Pak-McK criterion. When fitting the data 
with this expression, one obtains M cr st = 2.46 and 0.48 
for the simulations and the experiments, respectively. 
The comparison of the results from experiments, sim- 
ulations and the scaling argument thus show that Wi c 
at the onset of elastic instability in a serpentine channel 
scales as the square root of R/W with a small offset when 
R/W tends towards zero. 

We have also tested the nature of the instability in 



the experiments. First we have shown that increasing 
or decreasing Q does not lead to different values of Wi c 
and no hysteresis is observed. Second, we have perturbed 
the flow at the inlet using a modified serpentine channel 
composed of a series of loops of small radius 50 jim < 
Ri < 250 um, followed by a serpentine channel of larger 
radius, R2 = 450 um (see Fig. [5]). No modification of the 
threshold in the second serpentine channel is observed 
when varying the radius of the first channel. The streak 
image shown in Fig. [5] was taken at a flow rate where the 
flow in the first channel is unstable whereas the flow in 
the second channel is stable. One observes a fast decay of 
the perturbation once the fluid enters the second channel. 
We thus conclude that the instability in the serpentine 
channel is supercritical in contrast to observations made 
in Couette or plate- plate devices [14 1. 

Using a combined experimental and numerical ap- 
proach we investigated the geometrical scaling of a 
purely-elastic instability in a serpentine channel. In this 
way we have shown that the critical Wcissenberg number 
at the onset of the instability scales as the square root 
of the inner radius divided by the width of the channel, 
R/W , with a small offset when R/W —> 0, in agreement 
with a simple theoretical analysis based on the dimen- 
sionless criterion of Pakdel-McKinley [H, [H . 
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